#核函数
import random
import time
t=time.time()
def kernel(a,b,kernel_type=None,sigma=None):
    if kernel_type==None:
        return numpy.dot(a,b.T)
    elif kernel_type=="rbf":
        if a.ndim==2 and b.ndim==2:
            m1, n1 = a.shape
            m2, n2 = b.shape
            a=numpy.tile(a,(m2,1,1))
            b=numpy.tile(b,(m1,1,1)).transpose((1,0,2))
            return numpy.exp(-numpy.sum((a-b)**2,axis=2)/(2*sigma**2))
        elif a.ndim==1 and b.ndim==1:
            return numpy.exp(-numpy.sum((a-b)**2)/(2*sigma**2))
def clipAlpha(alpha,H,L):

    if alpha > H:
        alpha = H
    if L > alpha:
        alpha = L
    return alpha
def selectJrand(i,m):
    j=i #we want to select any J not equal to i
    while (j==i):
        j = int(numpy.random.uniform(0,m))
    return j
from matplotlib import pyplot
import numpy
numpy.random.seed(100)
data=numpy.loadtxt("./testSet.txt",delimiter="\t")
# pyplot.scatter(data[:,0],data[:,1],c=data[:,-1],edgecolors="red",linewidths=1.5,alpha=0.7)
# pyplot.show()
label=data[:,-1]
feature=data[:,:-1]
m,n=feature.shape
#给定罚因子
C=100
#初始化alpha和b
b=0
alpha=numpy.zeros(m)
#查找精度
tol=1e-5
#最大迭代次数
max_iters=400
#计算初始误差
E=numpy.sum(label*alpha*kernel(feature,feature),axis=1)+b-label
#开始选择优化变量,检查是否满足KKT条件
iters=0
flag=True
c=0
while c<max_iters:
    li=[]
    for k in range(m):#扫描整个数据集
        #第一个变量的选择
        if ((label[k] * E[k] < -tol) and (alpha[k] < C)) or ((label[k] * E[k] > tol) and (alpha[k] > 0)):
            li.append(k)
    #找到所有违反KKT条件的样本
            #第二个变量的选择
    idx = numpy.arange(m)[(alpha > 0) & (alpha < C)]
    if not li:
        break
    if flag:
        for i in li:
            if alpha[i]>0 and alpha[i]<C:

                # j=numpy.argmax(numpy.abs(E[li]-E[i]))#优先在非边界样本总查找
                # j=li[j]
                # if E[j]==E[i]:
                #     if len(idx)==0:
                #         j = selectJrand(i, m)
                #     else:
                #         j = idx[numpy.argmax(numpy.abs(E[(alpha > 0) & (alpha < C)] - E[i]))]
                #         if E[i] == E[j]:
                #             j = selectJrand(i, m)
                idx = numpy.arange(m)[(alpha > 0) & (alpha < C)]
                if len(idx) == 0:
                    j = selectJrand(i, m)
                else:
                    j = idx[numpy.argmax(numpy.abs(E[(alpha > 0) & (alpha < C)] - E[i]))]
                    if E[i] == E[j]:
                        j = selectJrand(i, m)
                    #continue
                #计算新的alpha的值
                #j = selectJrand(i, m)
                alphaj_old = alpha[j].copy()
                alphai_old = alpha[i].copy()
                eta=kernel(feature[i],feature[i])+kernel(feature[j],feature[j])-2*kernel(feature[i],feature[j])
                #print((E[i]-E[j])/eta)
                if eta<=0:
                    continue
                alphaj_new=alphaj_old+label[j]*(E[i]-E[j])/eta

                if (label[i] != label[j]):
                    L = max(0, alphaj_old - alphai_old)
                    H = min(C, C + alphaj_old - alphai_old)
                else:
                    L = max(0, alphaj_old + alphai_old - C)
                    H = min(C, alphaj_old + alphai_old)
                if L==H:
                    continue
                alpha[j] = clipAlpha(alphaj_new, H, L)
                #alpha[j]=alphaj_new

                alpha[i]=alphai_old+label[j]*label[i] * (alphaj_old - alpha[j])
                b1=b-E[i]-label[i]*kernel(feature[i],feature[i])*(alpha[i]-alphai_old)-label[j]*kernel(feature[i],feature[j])*(alpha[j]-alphaj_old)
                b2=b-E[j]-label[i]*kernel(feature[i],feature[j])*(alpha[i]-alphai_old)-label[j]*kernel(feature[j],feature[j])*(alpha[j]-alphaj_old)
                if (0 < alpha[i]) and (C > alpha[i]):
                    b=b1
                elif (0 < alpha[j]) and (C > alpha[j]):
                    b=b2
                else:
                    b=(b1+b2)/2
                E = numpy.sum(label * alpha * kernel(feature, feature), axis=1) + b - label
                flag=False
                print(b)
                c=c+1
                #break
        flag=False
    else:
        for i in li:
            # j = numpy.argmax(numpy.abs(E[li] - E[i]))
            # j = li[j]
            # if E[i]==E[j]:
            #     #j = selectJrand(i, m)
            #     if len(idx)==0:
            #         j = selectJrand(i, m)
            #     else:
            #         j=idx[numpy.argmax(numpy.abs(E[(alpha>0)&(alpha<C)] - E[i]))]
            #         if E[i]==E[j]:
            #             j = selectJrand(i, m)
            idx = numpy.arange(m)[(alpha > 0) & (alpha < C)]
            if len(idx) == 0:
                j = selectJrand(i, m)
            else:
                j = idx[numpy.argmax(numpy.abs(E[(alpha > 0) & (alpha < C)] - E[i]))]
                if E[i] == E[j]:
                    j = selectJrand(i, m)
                #continue
            # 计算新的alpha的值
            #j = selectJrand(i, m)
            alphaj_old = alpha[j].copy()
            alphai_old = alpha[i].copy()
            eta = kernel(feature[i], feature[i]) + kernel(feature[j], feature[j]) - 2 * kernel(feature[i], feature[j])
            # print((E[i]-E[j])/eta)
            if eta <= 0:
                continue
            alphaj_new = alphaj_old + label[j] * (E[i] - E[j]) / eta

            if (label[i] != label[j]):
                L = max(0, alphaj_old - alphai_old)
                H = min(C, C + alphaj_old - alphai_old)
            else:
                L = max(0, alphaj_old + alphai_old - C)
                H = min(C, alphaj_old + alphai_old)
            if L == H:
                continue
            alpha[j] = clipAlpha(alphaj_new, H, L)
            # alpha[j]=alphaj_new

            alpha[i] = alphai_old + label[j] * label[i] * (alphaj_old - alpha[j])
            b1 = b - E[i] - label[i] * kernel(feature[i], feature[i]) * (alpha[i] - alphai_old) - label[j] * kernel(
                feature[i], feature[j]) * (alpha[j] - alphaj_old)
            b2 = b - E[j] - label[i] * kernel(feature[i], feature[j]) * (alpha[i] - alphai_old) - label[j] * kernel(
                feature[j], feature[j]) * (alpha[j] - alphaj_old)
            if (0 < alpha[i]) and (C > alpha[i]):
                b = b1
            elif (0 < alpha[j]) and (C > alpha[j]):
                b = b2
            else:
                b = (b1 + b2) / 2
            E = numpy.sum(label * alpha * kernel(feature, feature), axis=1) + b - label
            #print(i,j)
            print(b)
            c=c+1
        flag = True
print(b)
print(time.time()-t)
print("OK")
pyplot.scatter(feature[:,0],feature[:,1],c=label)
pyplot.scatter(feature[:,0][alpha>0],feature[:,1][alpha>0],c="none",edgecolors="red",linewidths=1.5,alpha=0.7,s=150)
pyplot.show()
#优化的SMO算法